\documentclass{revtex4}
\usepackage{color}
\usepackage{colordvi}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{psfrag}
\usepackage{epsfig}
\usepackage{bm}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\rv}{\mathbf{r}}
\newcommand{\Rv}{\mathbf{R}}
\newcommand{\Lv}{\mathbf{L}}
\newcommand{\kv}{\mathbf{k}}
\newcommand{\qv}{\mathbf{q}}
\newcommand{\hv}{\mathbf{K}}
\newcommand{\Gv}{\mathbf{G}}
\newcommand{\Tr}{\mathrm{Tr}}
\newcommand{\Op}{\mathcal{O}}
\newcommand{\dtau}{\tau}
\newcommand{\imt}{t} %imaginary time
\newcommand{\barN}{\bar{N}}
\newcommand{\barS}{\overline{S}}
\newcommand{\barV}{\overline{V}}
\newcommand{\shah}{\amalg\!\!\amalg}

\begin{document}
\title{Specific heat in \sc{quest}}

\author{Simone Chiesa}%\email{}
\affiliation{W \& M}

\maketitle

We are interested in computing the specific heat of 
\begin{equation}
H=\sum_{ij\sigma} t_{ij}^\sigma c_{i\sigma}^\dagger c_{j\sigma}
+\sum_i U_i c^\dagger_{i\uparrow} c_{i\uparrow} c^\dagger_{i\downarrow} c_{i\downarrow}
\end{equation}
using $\langle H^2 \rangle - \langle H\rangle^2$. We break 
$ H^2 $ in four pieces $H^2 = TT + TU + UT + UU$ with
\begin{equation}
TT = \sum_{ij\sigma}\sum_{kh\alpha}  t_{ij}^\sigma t_{kh}^\alpha c_{i\sigma}^\dagger c_{j\sigma} c_{k\alpha}^\dagger c_{h\alpha}
\end{equation}
\begin{equation}
TU = \sum_{ij\sigma}\sum_{k}  t_{ij}^\sigma U_{k} c_{i\sigma}^\dagger c_{j\sigma} c^\dagger_{k\uparrow} c_{k\uparrow} c^\dagger_{k\downarrow} c_{k\downarrow}
\end{equation}
\begin{equation}
UT = \sum_{ij\sigma}\sum_{k}  t_{ij}^\sigma U_{k} c^\dagger_{k\uparrow} c_{k\uparrow} c^\dagger_{k\downarrow} c_{k\downarrow}  c_{i\sigma}^\dagger c_{j\sigma}
\end{equation}
\begin{equation}
UU = \sum_{ik} U_i U_k c^\dagger_{i\uparrow} c_{i\uparrow} c^\dagger_{i\downarrow} c_{i\downarrow} c^\dagger_{k\uparrow} c_{k\uparrow} c^\dagger_{k\downarrow} c_{k\downarrow}.
\end{equation}
For a given configuration of auxiliary fields, the expectation value can be computed via Wick's theorem.
To do so, it is useful to introduce the following quantities
\begin{eqnarray}
G_{ji}^\sigma  = \langle c_{j\sigma} c^\dagger_{i\sigma} \rangle \\
n_{i\sigma}  = \langle c^\dagger_{i\sigma} c_{i\sigma} \rangle \\
A_{ji}^\sigma  = \sum_k G_{ik}^\sigma t_{kj}^\sigma 
\end{eqnarray}
We then have (the notation closely mirrors the one used in the code with the exception of the matrix
$t$ which is called $h$ inside {\sc quest})
\begin{equation}
\begin{split}
\langle TT \rangle 
&=\sum_{ij\sigma}\sum_{kh\alpha}  t_{ij}^\sigma t_{kh}^\alpha \Big[
\langle c_{i\sigma}^\dagger c_{j\sigma}\rangle\langle c_{k\alpha}^\dagger c_{h\alpha} \rangle +
\langle c_{i\sigma}^\dagger c_{h\alpha}\rangle \langle c_{j\sigma} c_{k\alpha}^\dagger  \rangle \Big]
\\
&=\Big[\sum_{j\sigma} (t_{jj}^\sigma - A_{jj}^\sigma)\Big]^2 +
\sum_{hj\sigma} (t_{hj}^\sigma - A_{hj}^\sigma) A_{jh}^\sigma \\
\end{split}
\end{equation}

\begin{equation}
\begin{split}
\langle TU \rangle 
&=\sum_{ij\sigma}\sum_k t_{ij}^\sigma U_k 
\Big[\langle c_{i\sigma}^\dagger c_{j\sigma}\rangle 
\langle c^\dagger_{k\uparrow} c_{k\uparrow}\rangle
\langle c^\dagger_{k\downarrow} c_{k\downarrow} \rangle
+
\langle c_{i\sigma}^\dagger c_{k\sigma}\rangle 
\langle c_{j\sigma} c^\dagger_{k\sigma}\rangle
\langle c^\dagger_{k -\sigma} c_{k -\sigma} \rangle\Big] \\
&=\Big[\sum_k U_k n_{k\uparrow} n_{k\downarrow}\Big]
\Big[\sum_{j\sigma} t_{jj}^\sigma - A_{jj}^\sigma\Big]
+
\sum_{jk\sigma} U_k n_{k-\sigma} (t_{kj}^\sigma- A_{kj}^\sigma) G_{jk}^\sigma
\end{split}
\end{equation}

\begin{equation}
\begin{split}
\langle UT \rangle 
&=\sum_{ij\sigma}\sum_k t_{ij}^\sigma U_k 
\Big[\langle c_{i\sigma}^\dagger c_{j\sigma}\rangle 
\langle c^\dagger_{k\uparrow} c_{k\uparrow}\rangle
\langle c^\dagger_{k\downarrow} c_{k\downarrow} \rangle
+
\langle c_{k\sigma}^\dagger c_{j\sigma}\rangle 
\langle c_{k\sigma} c^\dagger_{i\sigma}\rangle
\langle c^\dagger_{k -\sigma} c_{k -\sigma} \rangle\Big] \\
&=\Big[\sum_k U_k n_{k\uparrow} n_{k\downarrow}\Big]
\Big[\sum_{j\sigma} t_{jj}^\sigma - A_{jj}^\sigma\Big]
-
\sum_{jk\sigma} U_k n_{k-\sigma} A_{kj}^\sigma G_{jk}^\sigma
+\sum_{k\sigma} U_k n_{k-\sigma} A_{kk}^\sigma
\end{split}
\end{equation}

\begin{equation}
\begin{split}
\langle UU \rangle 
&=\sum_{ik} U_i U_k
\Big[
\langle c_{i\uparrow}^\dagger c_{i\uparrow}\rangle 
\langle c_{k\uparrow}^\dagger c_{k\uparrow}\rangle +
\langle c_{i\uparrow}^\dagger c_{k\uparrow}\rangle 
\langle c_{i\uparrow}^\dagger c_{k\uparrow}\rangle 
\Big]
\Big[
\langle c_{i\downarrow}^\dagger c_{i\downarrow}\rangle 
\langle c_{k\downarrow}^\dagger c_{k\downarrow}\rangle +
\langle c_{i\downarrow}^\dagger c_{k\downarrow}\rangle 
\langle c_{i\downarrow}^\dagger c_{k\downarrow}\rangle 
\Big] \\
&=\sum_{ik} U_i U_k
\Big[
n_{i\uparrow} n_{k\uparrow} +
(\delta_{ik}-G_{ki}^\uparrow)G_{ik}^\uparrow
\Big]
\Big[
n_{i\downarrow} n_{k\downarrow} +
(\delta_{ik}-G_{ki}^\downarrow)G_{ik}^\downarrow
\Big]
\end{split}
\end{equation}


\subsection*{What gets computed in {\sc quest}}

{\sc quest} works in the grand-canonical ensemble with a specific heat computed as
\begin{equation}
C_\mu = \frac{1}{N_s}\frac{\partial\langle H - \mu N \rangle}{\partial T} = \frac{\langle (H - \mu N)^2 \rangle-\langle (H - \mu N) \rangle^2}{N_s T^2}
\end{equation}
where $N_s$ is the number of unit cell in the cluster.
Note that this corresponds to
\begin{equation}
\frac{1}{N_s}\left.\frac{\partial S}{\partial T}\right|_{\mu} = \frac{C_\mu}{T}
\end{equation}
{\em i.e.} $C_\mu$ is related to the derivative of $S$ at constant $\mu$, not constant density.
One can get to the latter using standard thermodynamic manipulations
\begin{equation}
\left.\frac{d S}{dT}\right|_{N} - \left.\frac{d S}{dT}\right|_{\mu} = 
-\left.\frac{\partial N}{\partial T}\right|_{\mu}^2
\left.\frac{\partial N}{\partial \mu}\right|_{T}^{-1}
\end{equation}
We are often interested in cases where $\partial N /\partial T =0$ ({\em i.e.} a particle-hole
symmetric Hamiltonian) in which case $C_\mu = C_N$.


\end{document}
